o 

(N 



(N 



O 

o 



X 



Time Evolution of the Spread of Diseases with a General 
Infectivity Profile on a Complex Dynamic Network 

Bahman Davoudi^ and Babak Pourbohloul^'^'* 

^Division of Mathematical Modeling, University of British Columbia Centre for 

Disease Control, Vancouver, British Columbia, Canada 

^School of Population & Public Health, University of British Columbia, Vancouver, 



May 25, 2010 



^ ■ British Columbia, Canada 

o ' * 

CD ; Corresponding author 

-)— > 

a 

B 

I . Abstract 

Pi I This manuscript introduces a new analytical approach for studying the time evolution of dis- 

ease spread on a finite size network. Our methodology can accommodate any disease with a 
general infectivity profile. This new approach is able to incorporate the impact of a general in- 
tervention - at the population level - in a number of different ways. Below, we discuss the details 

^ ' of the equations involved and compare the outcomes of analytical calculation against simulation 

^ . results. We conclude with a discussion of possible extensions of this methodology. 

O 

ir] 1 Introduction 

o 
o 

The time evolution of disease spread within human populations is a very interesting and multifaceted 
topic; it is closely related to the rate of new infections in a population, which has significant impli- 
cations when designing public health policy. The spread of disease within a population is complex 
and depends on a number of different factors, including, social connectivity patterns, cultural prac- 
tices and education surrounding hygiene and intervention strategies, the level of preexisting immunity 
and finally, the infectivity of the infectious agent. In order to produce a reliable estimate of the new 
infection rate, a model should incorporate - at the very least - the abovementioned factors. 

The spread of disease in a population involves a complex stochastic branching process. This pro- 
cess is comprised of three distinct constituents, namely, the stochastic phase, exponential phase and 
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declining phase. To date, each phase has undergone a considerable amount of scrutiny using two main 
techniques, i.e., compartmental and network models. The stochastic phase was studied using discrete 
and continuous time approaches [|Il|2l[3l|4l, while the exponential and declining phases were studied 
using a variety of models and techniques, including, compartmental JSl [6l |71 and network models 
[|2l |3l m [81 m. Both of these method, however, have their downfalls; that is, compartmental models 
deal solely with constant infection and removal rates, do not incorporate any memory of infections, 
distribute the force of infection uniformly in the population, and include the finite size eff"ect in a very 
simple manner. Moreover, some of the limitations associated with network methods are seen in the 
work by Nel et.al [(2l and Davoudi et.al [9\. They considered disease transmission to follow the gen- 
eration time concept, that is, they assumed infection or recovery/removal occurred at a discrete time 
within the period of r. This led to a significant simplification of the calculations, thus, the model be- 
came unrealistic for a disease with a long period of infection. The work by Volz et.al flSl included the 
finite size efl'ect in a comprehensive manner, but the method could not accommodate diseases with 
a complex infection profile. To circumvent these downfalls we present a new methodology which 
can sustain complex infection and recovery/removal profiles, and can take the finite size effect into 
account, in a precise manner. 

Below, we first discuss the required theoretical components for this methodology. We then com- 
pare theoretical results against simulation results to understand the precision of the method and the 
level of approximation involved. Subsequently, we discuss the possibility of extending this current 
methodology to address systems with a complex dynamic network structure. 

2 Theory 

2.1 Basic Notions 

We consider the preexisting contact network to be the platform for disease transmission within a 
population. Individuals are represented by a vertex and contact between two individuals is represented 
by a link. Links are composed of two stubs, each of which is attached to its own vertex. For a 
population of size A'^, we create a random network where the degree of connectivity between vertices 
is coded by the probability distribution function, p^, where k is the degree. The average excess 
degree lEl is defined as the degree of secondarily infected vertices and is given by Z^ = ziizi where 



Zi = (k)p, is the average degree, zi = {k^)p, - {k)p, and finally {k")p^ = 2^=0 ^"Pk- 

The process of disease transmission within a population is a complex phenomenon. Individuals, 
or vertices, are infected and removed at specific times, denoted by tt and t,- respectively. Infected 
vertices are removed from the population for a number of reasons, including death, quarantine, and 
recovery for example. Moreover, the age of infection and the time between infection and removal, or 
the "time-to-removal", for each infected vertex are defined as t = t - ti and t,. = t^ - U, respectively. 
The transmission of infection between an infected vertex and a susceptible vertex is dictated by the 
infectivity function, /1;(t), in which /l,(r)JT denotes the probability of transmission during the interval 
T and T -I- dr. In the same manner, we define the removal function as Ar{T), in which Ar{T)dT yields the 
chance of removal of an infected vertex during the interval r and r -i- dr. The transmissibility, T{t, t,.), 
provides the probability of transmission until time r, for an infected vertex that is removed after time 
Tr and satisfies the following equation lfTOlfTTI 

1 - exp i- L Ai(u)du) T < T,- 

T{T,Tr) = \ ^ (1) 

1 - exp (- r '' Ai(u)du] otherwise 
We define *F(tJ [.10J as the probability that a vertex has a time-to-removal > r,., which is given by 

¥(T,-) = expj- r Ariu)du\, (2) 

subject to the condition W(oo) = 0. We also define the probability density function as lA(rr) = — '^^ 

J "CO ^ 

if/(u)du). The basic reproduction number is given by "^o = —T, where T is the ultimate 

transmissibility T = J i//(t)T(t,t < t^)[|3|]. 

Finally, J{t) denotes the rate of newly infected vertices at time t. At a given time t, a fraction, 

*P(t), of infected vertices J(t - t) with the age of infection t remain infectious. Therefore, the total 

number of infectious vertices, at a given time, can be calculated by 

Ni(t) = f J{t- T)'¥iT)dT (3) 

Jo 

and the number of removed and susceptible vertices are given by 

N,{t)= f J(t-T)[l-nr)]dT, (4) 

Jo 

and 

NAt)=N-Nr{t)-Ni{t). (5) 

respectively. 



2.2 Disease Transmission Dynamics on a Network 

In the following section we introduce and discuss the set of equations we use to find the rate of new 
infections J(t) as a function of time. This calculation becomes possible once we combine the network 
aspects (vertices connectivity) with the disease status of each vertex. To elaborate, we start with one 
infectious vertex with excess degree Zv, and assume that it was infected at time ti = 0. With this 
knowledge, we then calculate the number of new infection that arose from this infected vertex, at the 
later time t = t, using 

dT(T, T < Tr) 

J(T) = Y(T)— ^^- -Z,„ (6) 

UT 

where ^^^^'^^"•' is the contribution of each link to disease transmission between time r and r + dr, 
given that the vertex was not removed by time t; the resulting contribution of Zv link is given by 
'''^^dr^^'^ ^x- The equation above is then multiplied by *P(t) in order to take the chance of removal into 
account. The total number of infections caused by the first infected vertex is given by TZ^ where 
T = JJ' ^'{T)^^^^^j^dT = C° iff(T)T(T,T < T,.)dT is the ultimate transmissibility, which yields the 
probability of infection along a link. 

Equation ^ can be easily extended to the initial phase of an epidemic, assuming that the excess 
degree of all vertices is the same. In general, the renewal equation for J(t) is as follows lfT2l [3]| 

Jit - T)^iT)^^ '-IZ.it, T)dT (7) 

dT 

The right hand side of the above equation gives the total number of transmitting links A{t) at time t, 
which leads to 7(f) « A{t) infections |l9l- Equation ([7]), where Z^it, r) ^t Z^ = zilzi, can be used before 
the finite size effect becomes important, which is a valid assumption while Nr{t) + Ni{t) «; Ns{t). In 
the limit Nr{t) + Ni{t) ~ Ns{t) , only a fraction of Z^ is used to connect the infected and susceptible 
vertices. An appropriate approximation for Z^it, t) is given below. 

Zx{t, t) is calculated in two steps. First, the typical degree of infected vertices is calculated during 
the process of disease transmission. Second, an estimate of the average number of links an infected 
vertex - with an infectious period of r - could have with susceptible vertices, at time t, is made. 

The first step is easily preformed for a random network [|2l ISl- A vertex is randomly picked 
and assigned to a collected class, while the probability that the chosen vertex has the degree k is 
given by <?i(l) = Pk- The function argument shows the number of collected vertices. The degree 
distribution of the uncollected vertices is given by pui^) = Pk- The expected degree of the first 



collected vertex is z(l) = zi and the average degree of the uncollected vertices is z(l) = zi- A 
second vertex is then chosen by picking a random stub. The probability that the second vertex has 
degree k is given by qkC^) = kpk/z(l), the degree distribution of the collected vertices is calculated as 
follows pk{2) = (qki^) + qk{2))/2 and the degree distribution of the uncollected vertices is specified 
by Pk(2.) = (Npk - 2pk{2))/{N - 2). The probability that the nth chosen vertex has degree k is qk{n) = 
kpk{n - y)lz{n - 1) and, in the same manner, the degree distributions of collected and uncollected 
vertices are given by 

n 

Pk(n) = ^ qkin)/n, (8) 

1=1 

Pk(n) = (Npk - npk(n))/(N - n), (9) 

where z(n) = YjZo ^Pk{n) is the average degree of uncollected vertices after n collections. We define 
z{n) = 2^0 kpk{n) as the average degree of collected vertices after n collections. The latter equations 
take the following form in the continuous limit 

dpk{n) Pk(n) ( k A QQN 



dn n \z{n) 



dn N - n\ z{n) 

The average degree of infected vertices that are infected between time t and t + dti% given by 
^^^ {NAt) + Ni(t) + J(t)dt]z(Nr(t) + Ni(t) + J{t)dt) - {Nr{t) + Nm~z{Nr{t) + Nm 

'^^'^ = 7m ^'^^ 



= z(Nr(t) + Ni(t)) + [Nr(t) + Ni(t)] "^^^"^ 



dn 



(13) 



n=N,{t)+N,{t) 



In the same manner, we find that the average degree of removed Zr(t), infectious Zi(t), and susceptible 
classes Zs(t) with the following formulas, respectively 

Zr(t) = z(Nr(t)), (14) 

,^, [Nr(t) + Nmmrjt) + Nm - NrityzjNrit)) 

'^^'^ = W) ' ^''^ 

Zs-(t) = z(Nm, (16) 

Z^t, r) can now be easily estimated. The number of stubs belonging to infected vertices, between 

time t-T and t-T + dr, is given by Zj{t - T)J{t - T)dT. The probability of one of these stubs connecting 

to the stub of a susceptible vertex, at time t, is given by Zs{t)N s{t) I Z\N , and thus 

r^ ,^ . ,, Mt)Ns{t) 

Z^{t, t) ~ Zj{t - t) -^. (17) 

ziN 



We also considered other approximation with different renewal equations 

dT{T, T < Tr) , 



Jit) 



f 

Jo 



J(t-T)^(T)- 



dT 



-Z,NAt)/NdT, 



(18) 



This new methodology provides a solution for SIR compartmental models, whereby the force of 
infection, ;S(t), and recovery rate, yir), are functions of r [4J. This is equivalent to the approximation 

3 Numerical results for the most simplistic network 

In this section we present the numerical result for three stylized networks, namely binomial pt = 
(TjP''(^ - P)^~* with p = zi/{N - 1) (zi = 10), exponential pt = {I - exp(-l/j^)) Qxp{-k/K) (k = 10) 
and bimodal, which are depicted in figure [B The size of the networks were A^ = 10, 000. To obtain 
a numerical solution with the current methodology, we first created two sequences, namely z{n) and 
z{n), for a specific network using equation dS]) and ^. We then recursively used the renewal equation 
^ to obtain the number of infections at a later time, meanwhile Zjit), Zs{t) and Zv(?, r) were calculated. 
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Figure 1: Three different degree distributions, binomial zi = 10, exponential a: = 10 and bimodal 
network. 



In figure [2l we compare the number of removed vertices for the new calculation (Analytical-N) 
and the compartmental model (Analytical-C), against simulations for a binomial network. Top left 
panel zi = 10, i, = 0.35 and Ar = 2; top right panel Ai = 0.2 and Ar = I; bottom left panel Aj = 0.5 and 



A,. = 2; and bottom right panel Ai = 0.3 and A^ = \. Both analytical approaches performed well for the 
binomial network, within this range of parameter values. The current approach slightly overestimated 
the final size for small 'Rq values and slightly underestimated the final size for large "^o values. 
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Figure 2: zi = 10, Aj = 0.25 and A,. = 2 left top panel. At = 0.2 and A,- = 1 right top panel. At = 0.5 
and Ar = 2 left bottom panel, and finally Aj = 0.3 and A,- = 1 right bottom panel. 

In figure [3l we compare the number of removed vertices for the new calculation (Analytical-N) 
and compartmental model (Analytical-C), against simulations for exponential network. Top left panel 
K = 10, Aj = 0.35 and A^ = 2; top right panel Ai = 0.2 and A,- = I; bottom left panel Aj = 0.5 and 
Ar = 2; and, bottom right panel Aj = 0.3 and Ar = I. The compartmental model was inferior to 
thebinomial network; for exponential networks, as we expect that Zj,(t, r) ~ Z^ will become a poor 
predictor for any network with a very wide degree distribution. 

Finally in figurelH we compare the number of removed vertices for the new calculation (Analytical- 
N) and compartmental model (Analytical-C), versus simulations for bimodal network. Top left panel 
Aj = 0.25 and Ar = 2; top right panel Aj = 0.2 and Ar = 1; bottom left panel Aj = 0.5 and Ar = 2; and, 
bottom right panel Aj = 0.3 and Ar = \. 

The small deviation for the current formalism and simulation for low and high Kq values is due to 
overestimating/underestimating the number of infections from the number of transmitting links. 
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Figure 3: k = 10, Ai = 0.25 and A,- = 2 left top panel, /I, = 0.2 and A,- = 1 right top panel, At = 0.5 and 
Ar = 2 left bottom panel, and finally Ai = 0.3 and Ar = I right bottom panel. 

The compartmental model fails to correctly capture the dynamics of the epidemics for populations 
in which a sizable portion of individuals have low degree, which turns to be an important characteristic 
of realistic human contact networks. 

4 Extensions 



In this section, we discuss different extensions of the methodology to address other models such as, 
open system models (in which the number of vertices is a function of time), dynamic network models, 
susceptible-infectious-removed-susceptible SIRS models, and multi-type systems (including systems 
with different age groups, each of which having different connectivity, infectivity, susceptibility, etc). 

4.1 Multi-type system 

Here we consider a system consisting of more than one type of nodes,where each node type has 
different inter- and intra- degree distributions, transmissibilities, and removal distributions. We define 
a set of node types and index them with superscript a. We first define A^^", A^ir), /if (r) and p" as 
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Figure 4: Ai = 0.25 and A,- = 2 left top panel, At = 0.2 and Ar = I right top panel, /I, = 0.5 and Ar = 2 
left bottom panel, and finally /l, = 0.3 and Ar = I right bottom panel. 

the number of vertices , new removal function, infectivity function and degree distribution of type a, 
respectively. We also define Nait) as the number of vertices from type yS in the classes a. Here we 
assume the connectivity between vertices, in different node types, is completely random. The current 
method can be further extended if there is any preference for intra- or inter-connectivity between 
vertices in different types. 

Our renewal equationdV]) takes the following trivial form 



J"(t) = J]fAt-T)^'^r) 



dTP"{r, T < Tr) 
dT~ 



Z^^"(t,T)di 



(19) 



Z^", however, is now a more complex function and must be calculated using the following steps; 
first, the collected and uncollected degree distribution for each type must be found. To do so we 
randomly pick a vertex. The probability that the first chosen vertex is in type a is n"(l) = N" /N, as a 
result the probability that the chosen vertex has degree k and comes from type a is given by <?^(1) = 
n"(l)p1. The degree distributions of collected and uncollected classes are given by p^(l) = ^^(1) and 
/?^(1) = (N"p'^-n"(l)p'^(l))/(N-n"(l)). Now the second vertex can be randomly chosen by selecting 
a random stub. The probability that the second vertex has degree k and belongs to type a is given by 



q"^(2) = n"{2)kpl{2)lz"{l) where rf{2) = {N" - n"{l))/{N - 1) and z"{l) = Zkkp"{l). The degree 
distribution of collected and uncollected classes are given by p'li2) = (<?^(1) + <?^(2))/(n"(l) + n"{2)) 
and p^(2) = (N"p"^ - {n"(l) + Tf{2))pl{2))l{N -n"{\)- n'\2)). It is thus plausible that the nth chosen 
vertex is in type a and has degree k with the probability q1{n) = n"(n)kp'^(n - l)/z"(n - 1). In the 
same manner 

Therefore, z"{t), z"{t), z"{t), and z1{t) are given by 

dr{n) 



zy) = rcA^rco + A^r (0) + [a^. (o + A^r (oi 



dn 



n=N?(t)+N'^{t) 



(22) 



z?(0 = r(N^(t)), (23) 
[A^;(0 + A^f (0]r(A^,'^(0 + Nf(t)) - N'^itTziN'^it)) 

'"^'^ = Nf^) ' ^2'^ 

z«(0 = z''(N:(t)), (25) 

where Nf{t) = f^ r{t - T)'i"'iT)dT, N^it) = f^ r{t - r) [1 - *I'"'(t)] Jt and N"^{t) =N- Nf{t) - N^{t). 

Finally, we have 

Zf(t,T)^zy-T)^^^^. (26) 

■^ Z[N 

As an example, we study an exponential degree distribution with k = 10, Aj = 0.2 and A^ = I. We 
use the multi-type framework to show how easily one can replicate the results obtained earlier for the 

similar degree distribution in section} } the current formalism. We devide the network to groups of 

vertices, each of which has a specifc degree where {a, /3} = {k, k'}. Our multi-type renewal equations 
is given by 

A crude approximation for the contact matrix, Z,x, is described below. The probability of a stub from 
a node type a connecting to a stub from node type jS is given by N^^z^^jNzi ■ This implies that the total 
number of links going from type a to j3 can be approximated by N"z"N^Zi/Nzi and consequently the 
number links per vertex is given by zIN^z^lNzi . 

In figure [5] we compare the result of the current calculation (Analytical-N) against one type 
(Analytical-Cl), multi-type (Analytical-Cn) and simulation models for the above-mentioned network, 
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whereby /I, = 0.2 and A^ = I. The excellent agreement between the two methods demonstrates that 
a network with a general degree distribution can be examined as a set multi-type system, within the 
current approximation. Both the multi-type framework and the current formalism, as discussed in 
sectior(3l yield similar levels of error when predicting the epidemic curve. Thus, one can use ei- 
ther approach; however, the multitype approach may become very expensive computationally for a 
network with a very wide degree distribution. 
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Figure 5: a: = 10, /I, = 0.2 and A^ = I. 



4.2 Open system 

In an open system, the number of vertices is a function of time. The previous set of equations still 
hold for an open system, however, we must now keep track of entering and exiting vertices in each 
class, as well as the corresponding change in degree distribution. For example, the number of vertices 
in a susceptible class can be calculated from 

dN,(t) 



dt 



Yjy^inA^^ - ^outA^, t)). 



(28) 



where Yliniout,s{^, is the rate of entry/exit of susceptible vertices with degree k at time t. The degree 
distribution of collected and uncollected vertices can be calculated from 



dpkin, t) dpkin, t) dn liinjAK t) - Tiouuh-ik^ 



dt 



dn dt 



+ 



n(t) 



dpk(n, t) dpkin, t) dn liin^k, t) - YloutA^, 



dt 



dn dt 



+ 



N{t) - n(t) 



(29) 
(30) 
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where n,v, ,>(A:, t) = n,„ ,(fc, t) = n,„^(fc, t). The first term on the right hand side of both equations is 
the contribution of collecting vertices, the second term arises from vertices entering or exiting the 
network. The partial derivatives of both equations are calculated from (fTOl) and (fTT)) respectively 

dpk(n,t) pt(n(t),t)( k \ 



dt n{t) \z{n{t)) 

dpt(n,t) Pk(n(t),t) (^ k 



dt N(t) - n{t) \ z{n{t)) 

4.3 Dynamic network 

As another possible extension, we consider networks where the degree of each vertex is a function of 
time. Dynamic networks are a simple example of an open system with the constrains 

^ liinAk, t) = Yj n-',-(^' ^) (33) 

k k 

where a = {s, i, r} is an index for susceptible, infectious and removed vertices. Accordingly, the 
outflow of vertices with a specific degree from a given class should be replaced by the same number 
of vertices but with diff'erent degrees. This is a consequent of the fact that infection is instantaneous 
and that a removed vertex remains removed. Ilin/out,a(k, t) could have complex dynamics as long as 
the above constraints are satisfied. 

4.4 SIRS model 

For the SIRS model, we first need to introduce the probability function, which specifies the chance of 
re-infection over time, once the infected vertex has recovered. This variable is generally a function of 
time and can be measured with respect to any infection time reference. We define the susceptibility 
function, As{t), in which /1,(t)Jt gives the probability of an infected vertex becoming susceptible 
again in the interval r and T+dT. We define *Pff/3(T.v) as a probability function which give the probability 
of the time of movement from disease state a to yS, moreover 



Jr-Tr 




"^apiTs) = exp - Aap{u)du (34) 
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The number of susceptible and removed vertices is given by 

N,(t) = N- f Jit- T)dT + \ J{t- t)[1 - *P„(r)]JT, (35) 

Nr{t) 

Jo 

pi 

Jit - T)V¥rs{T) - T,>(T)]rfT. (36) 



I J(t-T)dT+ I , 

Jo Jo 

= f Jit-T)[l-^'ir(T)]dT- r 7(? - T)[l - ¥„(t)]Jt 

Jo Jo 

-f 

Jo 



The degree distribution of collected and uncollected vertices is calculated as follows 

dpkjn, t) _ dpkjn, t) dn _ Ilout,r{k, t) 
dt dn dt n{t) 



(37) 



dpt(n,t) dpt(n,t)dn n,„,,(fc,0 
dt dn dt N(t) - n{t) 

where Ilout,r{^, t) is the outflow of recovered vertices to susceptible classes and in the same manner 

^in,s{k, t) is inflow of susceptible vertices from the recovered class; this calculation involves defining 

the rate of outgoing recovered vertices as 

Joutit) = j^Jj(t- r)(l - Y„(t))Jt. (39) 

We also define the rate of outgoing degree of recovered vertices as 

Zoutit) = - f J{t- T)Zj{t - T)(l - T„(r))JT. (40) 

dt Jo 

The expected degree of outgoing recovered vertices is then given by 

kourit) = f4^, (41) 

■Joutyt) 

and as a result 

'n.outAk, t) = n,-„,,(fc, t) = 6 k,k„„(t) Joutit)- (42) 

In practice kgutit) is not an integer function since it gives the average degree of new susceptible ver- 
tices; thus, one must properly distribute new vertices around k„ut{t) to ensure that the average degree 
of the system remains constant. 

5 Conclusions and Discussion 

The novel methodology outlined above allows us to evaluate the time evolution of disease spread on 
a network. Our methodology is able to accommodate diseases with very general infectivity profiles. 
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Additionally, this methodology can manage multi-type networks, dynamical networks, and SIRS sys- 
tems . The precision of this methodology depends on the accuracy of the kernel of the renewal equation 
for the infection rate, which will be the subject of future investigations. 

6 Acknoledgment 

BP would like to acknowledge the support of the Canadian Institutes of Health Research (grant nos. 
MOP-81273, PPR-79231 and PTL-97126 [Team Leader grant (CanPan II)]) and the Michael Smith 
Foundation for Health Research (Senior Scholar Funds). BD was supported by these grants. 

References 

[1] M. Marder, Phys. Rev. E 75, 066103 (2007). 

[2] P-A. Noel, B. Davoudi, L. J. Dub, R. C. Brunham and B. Pourbohloul, Phys. Rev. E 79, 026101 
(2008). 

[3] B. Davoudi, J. C. Miller, R. Meza, L. A. Meyers, D. J. D. Earn, B. Pourbohloul, Submitted to 
Jour. Theo. Bio. 

[4] J. C. Miller, B. Davoudi, R. Meza, A. C. Slim, B. Pourbohloul, Jour. Math. Bio., 262, 107 
(2010). 

[5] N. T. J Bailey, The Mathematical Theory of Infection Disease and its applications (Hafner Press, 
New York, 1975). 

[6] R. M. Anderson and R. M. May, Infectious Disease of Humans (Oxford University Press, Ox- 
ford, 1991). 

[7] F Brauer, P. van den Driessche and J. Wu, Mathematical Epidemiology (Lecture Notes in Math- 
ematics, Mathematical Biosciences subseries 1945, Springer, 2008). 

[8] E. Volz. arXiv:0705.2092v l: E. Volz andL. A. Meyers, arXiv:070521vl. 

[9] B. Davoudi, F Brauer, B. Pourbohloul, submitted. 



14 



[10] D. R. Cox and D. Oakes, Analysis of survival data (Chapman & Hall 1984). 
[11] M. E. J. Newman, Phys. Rev. E 66, 016128 (2002). 
[12] A. J. Lotka, Ann. of Math. Stat. 10, 1 (1939). 



15 



6000 



^ 4000 - 



2000 




— Analytical-N 

— Analytical-C1 

— - Analytical-C2 

Simulation 



20 
t 



40 



0.5- 



— Analytical-N 

— Analytical-C 

— Simulation 



n I I I I I I I I I I I I I I I I I I I I I I I I I I I I I L 

10 20 30 

t 



